#ifndef ORTHANTSAMPLER_HPP
#define ORTHANTSAMPLER_HPP
//#include "Data/Data.hpp"
//#include "Gibbs.hpp"
template<class Kernel, class Resample>
class Sampler{
public:
	Sampler(Kernel *K,Resample *R,Orthant::Orthant *D,int M, double C){
		_p=(*D).Get_p();
		mat Xt(M,_p);//p le nombre de parametre a obtenir de Density
		W=new double[M];
		X=Xt;
		_K=K;
		_R=R;
		_Z=0;
		_C=C;
		_boo=0;
		_D=D;
		mat m(_p,1);
		m.fill(0);
		mat Sig(_p,_p);
		Sig.fill(1);
		Sig=diagmat(Sig);
		_G=new Distribution::Gaussian(_p,m,Sig);
		kk=-1;
		_thres=C*M;
		_M=M;
		Normalize();
		y=X.t();
	//	cout << "///" << _n << "///";
	}
	~Sampler(){
	//	delete[] W;
	}
	void Filter()
	{
		this->Init();
		int i=1;
		double p=0; 
		//for(int i=1;i<(_n);i++)
	//	i{
		double b=(*_D).Get_gamma();
	//	b+=(double)(1/(double)(_P-1));
		_temp=StepLength(b);
		Phiv.push_back(_temp);
		b+=_temp;
//		cout << "b: " << b << " p";
		_D->Set_gamma(b);
		y=X.t();
		
		_D->Weight(y,W);		
		
		double ess=Ess_W(W,1);
		cout << " ESS: " << ess;
		X=y.t();
		p=_D->Get_gamma();
		while(p<1){
	//		cout << i;
			cout << " step: " << i << "\n";
			this->Step(i);
			p=_D->Get_gamma();
		}
	//	Correction2();
		
	}
	void Init()
	{
		//pior
		X=_G->scrambled(_M);
	}
	void Step(int i)
	{
		//cout << "Ess:" << Ess << "\n";
				//cout << "\\";
//		if(Ess< _thres)
//		{
			mat v=growingvect(_M);
			double *w=new double[_M];
			for(int j=0;j<_M;j++){
				 w[j]=exp(W[j]);
			//	cout << w[j] << "\n";
			}
			(*_R)(&_M,w,v);
			delete[] w;
			//cout << v << "\n";
			
			Arangemat(X,v);
			Normalize();//poids fixé a 1
			//mat C=cov(X);
			//cout << C;
			//cout << "sum: "  << sum(X,0) << "\n";
			//cout << "before:" << X.row(1);
	//		for(int k=0;k<3;k++)
	//		{	
			
			//cout << "\nmean " << mean(X,0) << "\n";
			int kk=0;
			double acc=1;
		//	cout << "\nmean " << X(0,span::all) << "\n";
/////////////////////////////////////////////////////////////////////////////////////
			for(int i=0;i<1;i++)
			{
				_K->Move(&(X),_D->Get_gamma());
			}
			//cout << "\nmean " << X(0,span::all) << "\n";
			kk++;	
		/*	if(Acc[kk]<0.15)
			{
				_K->Set_Sigma(0.75);
			}else if(Acc[kk]>0.45){
			//	cout << "Hello\n" << Acc[Acc.size()];
				_K->Set_Sigma(1.5);
			}	
			*/
			
			//		}
			//cout << "after:" << X.row(1);

//		}
				
			y=X.t();
			//cout << X << "////";	
			double b=_D->Get_gamma();
		//	b+=(double)(1/(double)(_P-1));
			_temp=StepLength(b);
			Phiv.push_back(_temp);
			b+=_temp;
			//cout << "b: " << b << "p";
			_D->Set_gamma(b);

			//cout<< "before:" << W[0] ;
			//cout<< "after:" << W[0] ;
			//dans le cas de tempering peut etre determiner
			_D->Weight(y,W);//le p designera comment on avance 
		//	cout << "after:" << y.col(1);
			X=y.t();
			double ess=Ess_W(W,1);
			//cout << " ESS: " << ess;
			//cout << " ESS: " << Ess_W(W,1);
	}
	double Ess_W(double *w, int tf)//renormalisation???
	{
		double *foo2= new double[_M];
		int sum2=exp(w[0]);
		for(int i=1;i<_M;i++)
		{
			sum2+=exp(w[i]);
		}

		if(tf){
			_Z+=log(sum2)-log(_M);
		}
		return (double)sum2;

	}
	double StepLength(double phi)
	{
		double u=10-phi;
		double l=0;
		double thres=0.00001;
		double e=2;
		double alpha=0.01;
		double eta=0;	
		int Bool=1;
		while((e>thres)&(l<1-phi))
		{
			mat y=X.t();
		//	cout << y; 
		/*	for(int i=0;i<_M;i++)
			{
			     	cout << W[i] << "\n";
			}*/
			eta=Ess_W(_D->Eta_star(alpha,W,y,Bool),0)/_M;
			if(eta>_C)
			{
				l=alpha;
				alpha=(double)(alpha+u)/2;
			}else{
				u=alpha;
				alpha=(double)(alpha+l)/2;

			}	
			e=abs<double>(u-l);
		//	cout << "\\\\\\ "<< eta <<" \\\\\\\\\n"; 
			Bool=0;
		}
		return min2<double>(u,1-phi);
	}


	
	mat Get_theta(void){return X;}
	void Normalize(void){
		for(int i=0;i<_M;i++){
			W[i]=0;
		}
	}
	mat   Get_W(void){
		double sum=W[0];
		mat C(_M,1);
		for(int i=1;i<_M;i++)
		{
			double t=log_add(sum,W[i]);
			sum=t;
		}	
		for(int i=0;i<_M;i++)
		{
			C(i,0)=exp((W[i]-sum)+log(_M));
		}
		return C;
	}
	double Get_Z(void){
		return _Z;}
	mat Get_Phiv(void)
	{
		int n=Phiv.size();
		mat res(n,1);
		for(int i=0;i<n;i++)
		{
			res(i,0)=Phiv[i];
		}
		return res; 
	}
	mat Get_Acc(void)
	{
		int n=Acc.size();
		mat res(n,1);
		for(int i=0;i<n;i++)
		{
			res(i,0)=Acc[i];
		}
		return res; 
	}

private:
	mat X;//value particle
	double *W;//weigths
	double Ess;
	int _M;
	double _thres;
	double _Z;
	Kernel *_K;
	Resample *_R;
	Orthant::Orthant *_D;	
	mat y;
	double _temp;
	int _p;
	int kk;
	vector<double> Phiv;
	vector<double> Acc;
	Distribution::Gaussian *_G;
	int _boo;
	double _C;
};



#endif
